/*
 * Copyright 2003 Sandia Corporation.
 * Under the terms of Contract DE-AC04-94AL85000, there is a non-exclusive
 * license for use of this work by or on behalf of the
 * U.S. Government. Redistribution and use in source and binary forms, with
 * or without modification, are permitted provided that this Notice and any
 * statement of authorship are reproduced on all copies.
 */

#include "vtkEdgeSubdivisionCriterion.h"
#include "vtkDataSetAttributes.h"
#include "vtkMatrix4x4.h"
#include "vtkStreamingTessellator.h"

#include <cmath>

void vtkEdgeSubdivisionCriterion::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
}

vtkEdgeSubdivisionCriterion::vtkEdgeSubdivisionCriterion()
{
  this->FieldIds = new int[vtkStreamingTessellator::MaxFieldSize];
  this->FieldOffsets = new int[vtkStreamingTessellator::MaxFieldSize + 1];
  this->FieldOffsets[0] = 0;
  this->NumberOfFields = 0;
}

vtkEdgeSubdivisionCriterion::~vtkEdgeSubdivisionCriterion()
{
  delete[] this->FieldIds;
  delete[] this->FieldOffsets;
};

void vtkEdgeSubdivisionCriterion::ResetFieldList()
{
  this->NumberOfFields = 0;
}

int vtkEdgeSubdivisionCriterion::PassField(int sourceId, int sourceSize, vtkStreamingTessellator* t)
{
  if (sourceSize + this->FieldOffsets[this->NumberOfFields] > vtkStreamingTessellator::MaxFieldSize)
  {
    vtkErrorMacro(
      "PassField source size (" << sourceSize << ") was too large for vtkStreamingTessellator");
    return -1;
  }

  int off = this->GetOutputField(sourceId);
  if (off == -1)
  {
    this->FieldIds[this->NumberOfFields] = sourceId;
    off = this->FieldOffsets[this->NumberOfFields];
    t->SetFieldSize(-1, this->FieldOffsets[++this->NumberOfFields] = off + sourceSize);
    this->Modified();
  }
  else
  {
    off = this->FieldOffsets[off];
    vtkWarningMacro("Field " << sourceId << " is already being passed as offset " << off << ".");
  }

  return off;
}

bool vtkEdgeSubdivisionCriterion::DontPassField(int sourceId, vtkStreamingTessellator* t)
{
  int id = this->GetOutputField(sourceId);
  if (id == -1)
    return false;

  int sz = this->FieldOffsets[id + 1] - this->FieldOffsets[id];
  for (int i = id + 1; i < this->GetNumberOfFields(); ++i)
  {
    this->FieldIds[i - 1] = this->FieldIds[i];
    this->FieldOffsets[i] = this->FieldOffsets[i + 1] - sz;
  }
  t->SetFieldSize(-1, this->FieldOffsets[this->GetNumberOfFields()]);
  this->Modified();

  return true;
}

int vtkEdgeSubdivisionCriterion::GetOutputField(int sourceId) const
{
  for (int i = 0; i < this->NumberOfFields; ++i)
    if (this->FieldIds[i] == sourceId)
      return i;

  return -1;
}

bool vtkEdgeSubdivisionCriterion::ViewDependentEval(const double* p0, double* p1, double* real_p1,
  const double* p2, int, vtkMatrix4x4* Transform, const double* PixelSize,
  double AllowableChordError) const
{
  double real_p1t[4];
  double intr_p1t[4];

  Transform->MultiplyPoint(real_p1, real_p1t);
  Transform->MultiplyPoint(p1, intr_p1t);
  double eprod = fabs(AllowableChordError * real_p1t[3] * intr_p1t[3]);
  /*
     fprintf( stderr, "eprod=%g, compare to <%g,%g>\n", eprod,
     fabs(real_p1t[0]*intr_p1t[3]-intr_p1t[0]*real_p1t[3])/PixelSize[0],
     fabs(real_p1t[1]*intr_p1t[3]-intr_p1t[1]*real_p1t[3])/PixelSize[1] );
   */
  if ((real_p1t[0] > real_p1t[3]) || (real_p1t[0] < -real_p1t[3]) || (real_p1t[1] > real_p1t[3]) ||
    (real_p1t[1] < -real_p1t[3]))
  {
    double p0t[4];
    double p2t[4];
    for (int i = 0; i < 3; i++)
    {
      p0t[i] = p0[i];
      p2t[i] = p2[i];
    }
    p0t[3] = p2t[3] = 1.;
    Transform->MultiplyPoint(p0t, p0t);
    Transform->MultiplyPoint(p2t, p2t);
    int p0Code = 0, p2Code = 0;
#define ENDPOINT_CODE(code, pt)                                                                    \
  if (pt[0] > pt[3])                                                                               \
    code += 1;                                                                                     \
  else if (pt[0] < -pt[3])                                                                         \
    code += 2;                                                                                     \
  if (pt[1] > pt[3])                                                                               \
    code += 4;                                                                                     \
  else if (pt[1] < -pt[3])                                                                         \
    code += 8;

    ENDPOINT_CODE(p0Code, p0t);
    ENDPOINT_CODE(p2Code, p2t);
    if (p0Code & p2Code)
    {
      return false;
    }
  }
  if (fabs(real_p1t[0] * intr_p1t[3] - intr_p1t[0] * real_p1t[3]) / PixelSize[0] > eprod ||
    fabs(real_p1t[1] * intr_p1t[3] - intr_p1t[1] * real_p1t[3]) / PixelSize[1] > eprod)
  {
    // copy the properly interpolated point into the result
    for (int c = 0; c < 3; ++c)
      p1[c] = real_p1[c];
    return true; // need to subdivide
  }

  return false; // no need to subdivide
}

bool vtkEdgeSubdivisionCriterion::FixedFieldErrorEval(
  double* p1, double* real_pf, int field_start, int criteria, double* AllowableL2Error2) const
{
  int id = 0;
  double mag;
  while (criteria)
  {
    if (!(criteria & 1))
    {
      criteria >>= 1;
      ++id;
      continue;
    }

    mag = 0.;
    int fsz = this->FieldOffsets[id + 1] - this->FieldOffsets[id];
    for (int c = 0; c < fsz; ++c)
    {
      double tmp = real_pf[c + field_start] - p1[c + field_start];
      mag += tmp * tmp;
    }
    if (mag > AllowableL2Error2[id])
    {
      return true;
    }
    criteria >>= 1;
    ++id;
  }

  return false;
}
